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Abstract. Sensitivity analysis can be a powerful tool for design optimization, error analysis 
and data assimilation in simulating engineering systems. For computing sensitivity of time averaged 
quantities in system exhibiting chaotic oscillations, conventional sensitivity analysis methods (tangent 
linear and adjoint methods) can produce erroneous results. This paper presents a new method for 
computing sensitivity derivatives for systems that exhibit both periodic and aperiodic (chaotic) 
limit cycle oscillations. This method overcomes the growth of sensitivity that disables traditional 
sensitivity analysis methods. 
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1. Introduction. Adjoint based computational sensitivity analysis of flow fields 
has become a powerful technology in designing airplanes ^1, 2 . In an adjoint based 
analysis, sensitivity to small perturbations is traced back from the quantity of inter- 
est, such as lift or drag coefficient, to perturbations to the flow field. The adjoint 
equations govern the backward propagation of sensitivities, and are derived from the 
governing equations of the flow. The solution to the adjoint equations reveals the 
gradient of the quantity of interest with respect to perturbations to the flow fleld, 
e.g., from modifying the surface geometry. This gradient can then be used to drive 
a gradient based optimization scheme to maximize or minimize the quantity of inter- 
est. The adjoint method has been shown to be able to solve very high dimensional 
optimization problems typically found in aerodynamic shape design and optimal flow 
control applications [31 H] . 

In unsteady aeroelastic problems, adjoint based methods can be a powerful and 
transformative technology. Sensitivity computed from such analysis can be used to 
design non-conventional airplanes, rotorcrafts and micro air vehicles to maximize 
performance and minimize cost and risk [3 |6l [71 [8]. In computing sensitivities of 
time-averaged quantities in unsteady aeroelastic simulations, many quantities of in- 
terest are time-averaged quantities, such as aerodynamic forces and structural fa- 
tigue. Despite recent advances in algorithms for sensitivity analysis unsteady fluid 
flows [9l Uni [m [71 [121 [13l [14] , sensitivity analysis of long time averages has shown 
to yield erroneous results for chaotic dynamical systems [15]. As illustrated in [15] 
and in Section [3] of this paper, naively computed sensitivity in a chaotic dynamical 
system grows as 0(e^"^) as the averaging length T increases, where A is the maximal 
Lyapunov exponent and has a positive value. 

The challenge of computing sensitivity of long time averages in chaotic dynamical 
systems has attracted much interest. The ensemble adjoint approach [15 [ I16 j approx- 
imates the sensitivity gradient by combining many adjoint evaluations on randomly 
sampled trajectories. However, the required number of trajectories can be very large, 
making the method computational costly even for small dynamical systems like the 
Lorenz attractor. Method based on linearizing the Fokker-Planck equation pT] com- 
putes sensitivity of long time averages by solving a PDE in the phase space. However, 
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a small random noise perturbation must be added to the dynamical system, whose 
effect is difRcult to quantify. Methods based on the Fluctuation-Dissipation Theorem 
has proven successful on some dynamical systems [181. However, they require the equi- 
librium distribution of the dynamical system to have properties similar to that of the 
Gaussian distribution. The recently proposed method based on Lyapunov eigenvector 
decomposition [19] computes accurate sensitivity on a short trajectory. However, the 
computational cost is linear to the number of positive Lyapunov exponents, making 
it expensive when many positive Lyapunov exponents coexist. 

This paper present a new computational method based on solving the least squares 
problem of the linearized equations. The solution of the least squares problem reveals a 
pair of shadowing trajectories. The algorithm is easy to implement. The computation 
cost of this new method is independent of the number of positive Lyapunov exponents, 
making it potentially suitable for high dimensional chaos. 

The rest of this paper is organized as follows. Section [5] describes the conventional 
methods for sensitivity analysis, including the tangent linear and the adjoint methods. 
Section |3] describes the challenge of using traditional sensitivity analysis methods for 
chaotic oscillations. Section 2] outlines our new method for computing sensitivity of 
both periodic and chaotic oscillations. Sections [SI H] and [7] apply this new method to 
the van der Pol oscillator, the Lorenz attractor and an aero-elastic oscillation model. 
Section [8] concludes the paper. 

2. Conventional method for sensitivity computation. Conventional tan- 
gent linear method for sensitivity computation [20] can apply to a autonomous dy- 
namical system 

^it,s)^fiuit,s),s) (2.1) 

parameterized by s, which can represent control variables, design variables, or uncer- 
tain parameters. The initial condition of the dynamical system is given as 

u{0,s)=uo{s) . (2.2) 

Conventional tangent linear method [20] of sensitivity analysis computes the derivative 
of an objective function 

J(s) = / J{u{t,s),s)dt (2.3) 



to the parameter s. 

The conventional tangent linear method starts by solving an initial value problem 
of the tangent equation 

dv^y df Of duo 



where Viv represents the linearized state 

' " ds 



v,y{t,s) ^ — it,s) (2.5) 



^The conventional tangent and adjoint methods work for non-autonomous systems. However, we 
focus on autonomous systems in this paper. 
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of the initial value problem. The desired sensitivity derivative can then be computed 
via chain rule: 

dJ fdJ dJ\ , 

ds Jq \ou ds J 

This method is also known as the forward sensitivity analysis method. Another 
popular method for sensitivity computation is the adjoint method |20| . It is the 
formal dual of the tangent linear method. The adjoint method has the advantage of 
being able to simultaneously compute the derivative of an objective function to many 
parameters. For this reason, it is often used in optimal control, inverse problems and 
data assimilation. 

3. Breakdown of the conventional method for time averaged output 
quantities. For a dynamical system exhibiting periodic and chaotic limit cycle os- 
cillations, many output variables are long time averaged output quantities, which 
represent long term behavior of the dynamical system. 

(J)(,s)= hm i r Jiuit,s),s)dt, (3.1) 

T^OO T Jq 

Under the ergodic assumption, such infinite time averages are independent of the 
initial condition u(0,s). They only depend on the parameter s. These quantities 
are often called the "climate" of the dynamical system. Many output quantities of 
key interest in aeroelastic oscillations, such as mean magnitude of oscillation and 
structural fatigue, can be approximated by long time averages. 

In a computational simulation, an infinite time average must be approximated 
with a finite time average 

1 '•^ 



(J)(s) « Jis) - - J{u{t, s), s) dt . (3.2) 

Because J(s) is in the form of Equation ()2.3p . the conventional method described 
in Section [2] can be used to compute dJ{s)/ds. One might think that since J is a 
good approximation to (J) when T is large, the computed sensitivity derivative dJ /ds 
should be a good approximation to the desired sensitivity derivative d(J)/ds when T 
is large enough. However, this is not the case. 

The conventional sensitivity analysis method described in Section 2 breaks down 
for long time averaged output quantities because dJ(s)/ds can be a terrible approx- 
imation to d{J)(s)/ds. They can be orders of magnitude different for large T. The 
finite time average approximation contains a small error J(s) — (J)(s). The magni- 
tude of this error is small for large T; however, its derivative to s can increase as T 
increases. As a result, the derivative dJ{s)/ds does is not a useful approximation of 
d{J){s)/ds, rather it grows exponentially with increasing T |15j . If one decompose 
J{u(t, s), s) into its time average and a zero-mean fluctuating term 

J{u{t, s),s) = ( J)(s) + J(u(t, s),s), 

the sensitivity derivative of J(s) to s can be written as 

dl{s) _d{J){s) ^ 1 ['^ dJ{u{t,s),s) 



ds ds T ,/o ds 
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while the first term is the desired sensitivity derivative, the second term often does 
not decay to as T increases [TS] [TT]. In particular, for chaotic dynamical systems, it 
is observed that du{t, s)/ds grows exponentially as t increases, due to the "butterfly 
efl'ect"[21 , also known as "chaotic error growth" jl5|. As a result, 

dJ{s) d{J){s) _ I ['^ dJ{u{t,s),s) I dJ dJdu{t,s) 



ds ds T Jq ds T Jq ds du ds 

increases exponentially as T increases. The dJ{s)/ds computed with either the con- 
ventional tangent linear method or the adjoint method diverges ft'om the desired 
sensitivity d{J){s)/ds. 

This divergence behavior is examined by Lea et al. |15| in studying the Lorenz 
system 

^^aiy-x), ^^^ip-z)^y, |=xy-/3z. (3.3) 

They observe that when T increases, the difference between dz/dp and dz/dp does 
not decrease but increases exponentially fast. Following Lea et al.'s numerical ex- 
periment, we compute z by time averaging z computed from the initial condition 
X = l,y = l,z = 30 over a time period of T = 2.26 and 131.4, respectively. The 
sensitivity derivative of z is computed with the procedure in Section [21 where the 
discrete linearization is used. Figure [01 shows the sensitivity derivative of the finite 
time averages for different averaging lengths T. The divergence of sensitivities is a 
fundamental challenge in computing sensitivities for time averaged output quantities 
in unsteady oscillations. 

4. Least Squares Sensitivity Analysis of time averaged output quanti- 
ties. Our new method takes advantage of the ergodic assumption of the dynamical 
system. The linearized equation ()2.4|) approximates the difference between two ad- 
jacent trajectories, one with a perturbation in the right hand side, and one without. 
The time averaged objective function we try to capture does not depend on the initial 
condition. Therefore, the initial condition of the linearized equation, describing the 
initial difference between the two adjacent trajectories, can take an arbitrary value. 
Our new method solves the linearized equation (|2.4|) without an explicit initial condi- 
tion. The freed-up degrees of freedom resulting from the lack of initial condition are 
used to explicitly control the growth of the solution of the linearized equation, thus 
suppressing the butterfly effect. 

4.1. Computing sensitivity of time averaged output from tangent direc- 
tion of a shadovkfing trajectory. Consider a semi-infinite trajectory u{t), t e [0, oo) 
satisfying 

/(»,.) («) 

at a parameter s, where / is sufficiently smooth with respect to u and s. Consider an 
infinitesimal perturbation = s + e. The solution u at the unperturbed parameter 
value s is no longer a solution at the perturbed parameter value Sg, but rather an 
approximate solution. By the shadowing lemma [22 , one can find a u^(t),t G [0,oo) 
that satisfies 



(4.2) 
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(c) 2 for T = 131.4 
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(b) Sensitivity of 2 for T = 2.26 
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(d) Sensitivity of J*^' for T = 131.4 



Fig. 3.1: Plots created with similar procedure to Lea et al|!l5 (permission granted). 
Left column: finite time average z plotted against a system parameter r; the infinite 
time average z has a slope of about 1. Right column: the sensitivity of the finite time 
average with respect to r, computed with method outlined in Section 2. As the length 
of time averaging T increases, these sensitivities do not approximate the sensitivity 
of the infinite time averaged dz/dr^ which is about 1. Note the order of magnitude of 
the y-axcs. 



such that u and Ue are infinitesimally close under an infinitesimal time transformation 
t ^ t + ea{t). Mathematically, both 



\u,{t + e(j{t)) 



and 



da 
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are uniformly bounded for t G [0,oo). Therefore, we can define the tangent direction 
of the shadowing trajectory 

^/^^^ ^ u^jt + eajt)) -u{t) 

By taking its time derivative, substituting Equations (I4.ip and (|4.2I) . using chain rule, 
and ignoring O(e^) terms, we obtain the governing equation for u': 

du' df , df da ^ ^ , , 

where da /dt is the rate of time dilation. 

By enabling uniform convergence and exchange of limit and derivative jl9| , the 
shadowing trajectory enables computing derivative of long time averaged quantities 
(J) as in Equation Let the finite time average of the perturbed solution be 



1 



T+e(T(T) 



T + ea{T) Jq 

With change of variable t — t + ea{t), we obtain dr = {1 + e da/dt) dt, and 

J. ^ ^.^^ [ JM + eait)ls) (l + 4) dt . 

The difference between Je and the unperturbed average J = ^ J{u, s) dt can be 
obtained via subtraction, 



Je - J dJ , dJ da _^ da— „ , , ,^ ^, 

where the over-line denotes time averaging in [0,T]. 

It is worth noting that a uniformly shadowing trajectory u' exists for the whole 
t 6 [0, oo) only when the attractor of the dynamical system is uniformly hyperbolic 
[23] |24j . A semi-infinite shadowing trajectory may not exist for quasi- hyperbolic 
and non-hyperbolic systems. However, the least squares procedure described in the 
next subsection often computes the tangent direction of a shadowing trajectory for 
sufficiently long time to permit accurate computation of sensitivities. 

4.2. Least squares for computing the tangent direction of a shadow- 
ing trajectory. The tangent direction of the shadowing trajectory u' we solve the 
constrained linear least squares system 

/"^ 1 T "^^ 2 7 dv df df 

mm / -V V + —rj dt , s.t. — ^ —v + — + rif . (4.5) 



j)(t),D(t) Jo 2 2 ' ' ' dt du ds 

The solution v{t) and rjit) in this problem approximate the tangent direction of the 
shadowing trajectory u'{t) and the rate of time dilation da/dt, respectively, v ap- 
proximates u' because their difference satisfies the equation 

d , „ df , ( da\ 
— iv — u) + \ n — 
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By decomposing v ~u' into Lyapunov eigenvectors v — u' — ^ cik{t)4>k(t), where each 
Lyapunov eigenvector satisfies 

at ou 

where Afc is the corresponding Lyapunov exponent. In particular, the Lyapunov ex- 
ponent corresponding to the zero Lyapunov exponent A^q = is (j)^^ = f{u). By 
substituting the decomposition oi v — u' into the We find that the coefficients a/j of 
the decomposition satisfy the equation (j4.6|) . we obtain 



dak _ I Afcflfe Afc ^ 

dt ]T] — da/dt Afc = 

For negative Afc, at decays exponentiahy, thus ak{t) ^ afc(O) for t ;» l/|Afc|. For 
positive Afc, a^ decays exponentiahy as t goes backward from T, thus ak(t) <^ ak{T) 
for T - i > l/|Afc|. Therefore, v - u' is smaU for t e [Tb,T - Te] other than a time 
transformation corresponding to a shift in the direction of ^fc,, = f{u). The buffer 
times are 

Tb = TTT"; TTTT i Te = 



max{Afc|Afc < 0} ' min{Afc|Afe > 0} ' 

When T 3> Tb + 7b, the desired sensitivity derivative can be approximated by sub- 
stituting V for u' and rj for da / dt in Equation (j4.4p . i.e.. 



d{J){s) dJ dJ — - 

4.3. Numerical solution of the least squares problem. To solve the con- 
strained least squares system, we first convert it from a continuous-time problem to a 
discrete-time problem. Here we use the trapezoidal time integration rule to discretize 
d/dt. For a given parameter value s, suppose the solution of the nonlinear equation 
is obtained at u„ = u{nAt),n = 0, . . . ,N ~ T/At. Denote the discrete solution as 
Vn = v{nAt) and rjn — rj{{in — i)At). The continuous least squares problem (14. 5p can 
be discretized as 

N II 119 N 



min$:^-,.^^f , s.t. 

n=0 n=l (4.: 



V OL > — 

2 ^2 

n— 1 



At 



where 



An — 'q^ {Um s) 



+ bn+Vnfn , I <n< N 



1 f df df , . 

On = :r -1-(Un,s) + -^(Un-l,s) , (4.9) 



In 



2 \ds ds 

f{Un,s) + /(w„_i,s) 
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The optimality condition of this hnear constrained least-squares problem can be rep- 
resented by the foUowing KKT system 



" / 








a" 


fl 
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where the matrices K, and G„ are 



m 

VN 



Wl 
WN 



71 = 0, , 



-b2 



(4.10) 



(4.11) 



This KKT system can be solved by Gauss eUmination of the lower-left block, thereby 
forming the Schur complement of the KKT system 



BB^w = b 



where 



B 



Fq 



Gi 
Fl 



G2 



Fn-1 ^ G 



N 













W2 




b2 


, w = 




, b = 






_WN _ 




Bn _ 



(4.12) 



(4.13) 



This symmetric positive definite Schur complement matrix BB"'" is block-tri- 
diagonal, where the block size is the dimension of the dynamical system m. It can 
be solved with a banded direct solver, which takes 0{N m?) floating point operations 
[25] . One can also perform a sparse QR factorization on the block-bi-diagonal B-'", 
and then use backward and forward substitution to compute w. The sparse QR fac- 
torization also takes 0{N m?) floating point operations f25j. In addition, iterative 
solution methods can potentially be used to efhciently solve the system even for large 

TO. 

Once m is computed, w„ and rjn can be computed from the upper blocks of the 
KKT system (|4.10p . By setting wq = w^+i = 0, we obtain 



= F^Wn+1 + Gf^Wn 



0,...,N 



fnWn 



1, 



.,N (4.14) 



Once Vn and rjn are computed, the desired sensitivity derivative can then be computed 
by discretizing Equation (|4.7|) 



d{J) 
ds 



N 

E 



T 
9n 



1 dJ 



(m„ 



(4.15) 
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where 

9n 



1 dJ 
N + l'du 



(u„, s) , n = 0, . . . ,iV , 



_ J_ / JjUn, s) + J(Mn_i, s) 

" " iV 2 



1 \ 

— ^ J(u„s) 

4=0 / 



(4.16) 



n = 1, 



4.4. Summary of the algorithm. 

1. Choose a sufficiently small time step size At and sufficiently large time domain 
size NAt. 

2. Compute a solution to the equation (|4.1I) at u„ = u(nAi), n = 0, . . . , iV. 

3. Compute the vectors and matrices /„, b„, Fn and Gn as defined in Equations 
(gH) and (|4TI|) . 

4. Choose a nonzero a, form matrix B and solve system (|4.12|) for w. 

5. Compute w„ and rjn from Equation (j4.14l) . 

6. Compute desired sensitivity using Equation (j4.15p . 

The computational cost of algorithm is 0{N ra^), where is the number of time 
steps and m is the dimension of the dynamical system. 

4.5. Adjoint formulation of the sensitivity computation method. An 

alternative method of computing the same sensitivity estimate as in Equation (j4.15p 
is through the corresponding discrete adjoint formulation. In this formulation, one 
first solves the adjoint system 



/ 
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fl 
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Vn 
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WN 



9o 
hi 
91 
h2 

92 

h-N 

9N 



(4.17) 



where the sub-matrices Fn, Gn and vectors /„, 6„ are defined as in Equations (|4.11|) . 
(|4.9p and (|4.16p . This system (|4.17p has exactly the same matrix as the system (|4.10p , 
but with a different right hand side. It can be solved by inverting 



BB^w = Bg , 



(4.18) 



where B is defined in Equation (|4.13p . w = {lii , ... ,111 n), and g is the upper portion of 
the right hand side of Equation (|4.17p . Once w is computed, dJ/ds can be computed 
via 



dJ 

ds 



N 1 ^ dJ 

n— 1 n— 



(4.19) 
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where 6„ is defined in Equation (I4.9p . It can be shown through matrix manipulation 
that the adjoint sensitivity estimate is equal to the sensitivity estimate in Section 2] 
up to round-olF error. The examples in the following sections use the algorithm in 
Section 

5. Application to the Van der Pol oscillator. We apply our method to the 
Van der Pol oscillator 

— ^-y + Pil-y)-. (5.1) 

The parameter in the system is /3. Figure [S.lal shows the limit cycle attractor as /3 




°S.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 

,a 

(c) dJ/d/3 as a function of /3 computed with fi- 
nite difference with step size A/9 = 0.05. The 
black dots are computed on 20 trajectories of 
length 50 at each value of /3. The red line is 
computed on a trajectory of length 5000 for each 



°S.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 

,a 

(d) dj/dp as a function of /3 computed with 
Least Squares sensitivity analysis. The black 
dots are computed on 20 trajectories of length 
50 at each value of /3. The red line is computed 
on a trajectory of length 5000 at each /3. 



Fig. 5.1: Least Squares Sensitivity Analysis of the van der Pol oscillator. 



varies from 0.2 to 2.0. As /? increases, the maximum magnitude of dy/dt significantly 
increases. Our objective function is chosen as the norm of dy/dt, which has 
a similar trend to the L°° norm and reflects the magnitude of the peak in dy/dt. 
Denoting u — {u^-^\ u'^^) = {y, dy/dt) as the state vector, we convert the second order 
ODE (|5.ip into two coupled first order ODEs. The objective function can be written 
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as 

2_ 1 

The method described in Section |4] is then appUed to compute v. for each /3, 
we start the simulation at t ^ —50, with uniform [0, 1] random initial conditions for 
both variables at t = —50. The ODE is integrated to t = to ensure that u{0) is 
almost on the limit cycle attractor. A trajectory u{t),0 < i < 50 is then computed 
using a scipy[26j wrapper of Isoda^Tj, with time step size At — 0.02. The resulting 
+ 1 = 2501 states m„ along the trajectory are used to construct the linear system 
(j4.10p for computing u„. 

The linearized solution v is then used to estimate the sensitivity of the output J 
with respect to the parameter l3: 



dp 50 7o 



(5.3) 



The computed sensitivity is compared against finite difference values in Figure 15.11 
Both finite difference and least squares sensitivity analysis is repeated 20 times on ran- 
domly initiated short trajectories and a long trajectory for each value of /3. The error 
in these calculations can be estimated from the spread of derivative values computed 
on the same parameter value. Figure 15.11 indicates that the least squares method is 
more accurate than finite difference for the same trajectory length. 

6. Application to the Lorenz system. We apply our method to the Lorenz 
system (|3.3p . and analyze sensitivity to the parameter p in the system. The behavior 
of the Lorenz system as p changes from to 100 is shown in Figure IS.lal and can be 
summarized as following |28) : 

• Stable fixed point attractor at (0, 0, 0) for < p <= 1. 

• Two stable fixed point attractors aX x — y = ±^ I3{p — 1), z = p — 1 for 
1< p< 24.74. 

• Quasi-hyperbolic strange attractors for 24.06 < p < 31. This includes the 
classic Lorenz attractor at p = 28. 

• Non-hyperbolic quasi-attractors for 31 < p < 99.5. 

• Periodic limit cycle attractors with an infinite series of period doubling for 
p > 99.5. 

Despite the many transitions in the fundamental nature of the system, the mean z 
value 



(z) = lim - / zdt (6.1) 



T^'^ T Jo 

apparently increases as the parameter p increases, (z) is chosen to be our time aver- 
aged output quantity in this study. 

Denoting u — (x, y, z), the method described in Section [4.4l is applied to compute 
v: for each p, we start the simulation at t = —50, with random uniform [0, 1] as 
initial conditions for x, y and z. The Lorenz equation is integrated to t = to ensure 
that u(0) is almost on the attractor. A trajectory u(i),0 < t < 50 is then computed 
using a scipvpS] wrapper of lsoda[37], with time step size Ai = 0.01. The resulting 
5001 states along the trajectory are used to construct the linear system (|4.10p for 
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(a) Attractors of the Lorenz system for p = 10 (b) z as a function of r, computed on 20 trajec- 
(open circle), p = 25, 50, 75 and 100 (blue, green, tories of length 50 for each value of p. 
red and black lines, respectively) 













i 

















°'°0 20 40 60 80 100 



P 

(c) dJ /dp as a function of p computed with finite 
difference with step size Ap = 2. The black dots 
are computed on 20 trajectories of length 50 at 
each value of p. The red line is computed on a 
trajectory of length 5000 at each p. 



1.4 - 




^0.8 . ••• 

0.6 
0.4 - 
0.2 - 

°°0 20 40 60 80 100 

P 

(d) dj/dp as a function of p computed with 
Least Squares Sensitivity Analysis. The black 
dots are computed on 20 trajectories of length 
50 at each value of p. The red line is computed 
on a trajectory of length 5000 at each p. 



Fig. 6.1: Least Squares Sensitivity Analysis of the Lorenz system. 



computing v. The computed hnearized state v is plotted with the primal state u in 
Figure [611 for four different values of p. The different cases correspond to a fixed point 
attractor {p = 10), a quasi- hyperbolic strange attractor {p — 25), a non- hyperbolic 
attractor (p — 75) and a limit cycle attractor [p = 100). 

The linearized solution v is then used to estimate the sensitivity of the output 
(J): = (z) with respect to the parameter p: 

The computed sensitivity is compared against finite difference values in Figure 16.11 
The dip in the finite difference value at around p ~ 22.5 is due to the discontinuous 
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transition from fixed point attractors to strange attractors at 24.0 < p < 24.74 (note 
that within this small range, the two types of attractors co-exist). For 24.74 < p < 
31, the Lorenz system is dominated by a quasi-hyperbolic attractor. Least squares 
sensitivity analysis computes accurate and consistent gradients on randomly chosen 
short trajectories on the attractor. 

As p increases beyond 31, the system is non- hyperbolic and its trajectories form 
an object known as a quasi-attractor [29]. For p > 99.5, the system transitions to 
periodic oscillations, then goes through an infinite series of period doubling bifurca- 
tions. The shadowing lemma on which the least squares sensitivity analysis is based 
no longer holds in these situations. However, we find that the method can still com- 
pute derivatives more accurate than finite difference method on the same trajectory 
lengths. 

7. Application to an aero-elastic limit cycle oscillator. We apply our 
method to a simple model of aeroelastic limit cycle oscillation, as shown in Figure mi 
The model is described in detail by Zhao and Yang [30]. The governing equations are 

-h 




Fig. 7.1: Model aero-elastic oscillator 



(Ph da dh 

—rr + 0.25 — + 0.1 — + 0.2 + 0.1 a = 

dt-' dt dt (7 1) 

d^h d^a da q „ 

0.25 — ^ + 0.5 — ^ + 0.1 — 0.5 a -f 20 - 0.1 Q a = 
dt'^ dt^ dt 

where h is the plunging degree of freedom, and a is the pitching degree of freedom. 
We analyze sensitivity to the reduced dynamic pressure Q. 

The bifurcation diagram of a as Q increases from 8 to 16 is shown in Figure [7. 2al 
The behavior of the system as Q varies is complex [3T]: At low values of Q, the system 
exhibits asymmetric periodic limit cycle oscillations. As Q increases beyond about 
10.25, a series of period doubling bifurcations occurs, leading to transition into chaos 
just beyond Q — 11. At about Q — 12.5, the system transitions to symmetric periodic 
limit cycle oscillation. At above about Q — 13.25, there appears to be small windows 
in which the limit cycle oscillation transitions discontinuously back to asymmetric 
periodic shape. Finally, at about Q = 13.9, the system recovers symmetric periodic 
limit cycle oscillations. The phase plot of the system at several values of Q is shown 
in Figure [7.2bl These include an asymmetric periodic limit cycle oscillation at Q = 8, 
a chaotic limit cycle oscillation at Q = 12, and a symmetric periodic limit cycle 
oscillation at Q = 16. 

We observe that the magnitude of the oscillation grows as Q increases. Here we 
choose the norm of the pitch angle a as the objective function. The norm has 
similar trend as the L°° norm, which indicates the magnitudes of the oscillations in the 



14 



Q. WANG AND R. HU 





8 9 10 11 12 13 14 15 16 

Q 

(c) dJ/dQ as a function of Q computed with 
finite difference with step size AQ = 0.2. The 
black dots are computed on 20 trajectories of 
length 300 at each value of Q. 



8 9 10 11 12 13 14 15 16 

Q 

(d) dJ/dQ as a function of Q computed with 
Least Squares Sensitivity Analysis. The black 
dots are computed on 20 trajectories of length 
300 at each value of Q. The red line is computed 
on a trajectory of length 30000 at each Q. 



Fig. 7.2: Least Squares Sensitivity Analysis on the aero-elastic oscillator model (|7.1|) 



pitching degree of freedom. Denoting u — {u^^\u^^\u^'^\u'^'^'') = (jj,a,dy/dt,da/dt) 
as the state vector, we convert the pair of second order ODEs (|7.ip into a set of four 
first order ODEs. The objective function can be written as 

(J) = (^lim 1 £ dt^ ' . (^lim 1 £ .(^) « dt^ ^ (7.2) 

We use the method described in Section l44l to compute the sensitivity derivative of the 
objective function to the input parameter Q. For each Q, we start the simulation at 
t = —300, with uniform [0, 1] random initial conditions for all 4 variables at t = —300. 
The ODE is integrated to i = to ensure that u(G) is almost on the attractor. A 
trajectory u{t), < t < 300 is then computed using a scipy[26 wrapper of lsoda[27]. 
with time step size At — 0.02. The resulting 15001 states along the trajectory are 
used to construct the linear system (j4.10p for computing v. 

The linearized solution v is then used to estimate the sensitivity of the output J 
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with respect to the parameter Q: 



d{J) J 



-:-7 



300 



U 



(7.3) 



dQ 



300 



The computed sensitivity is compared against finite difference values in Figure 17.21 
The sensitivity derivative computed using least squares sensitivity analysis is more 
accurate than finite difference when the system exhibits periodic and chaotic limit 
cycle oscillations. 

8. Conclusion. We presented a computational method for sensitivity analysis 
of long time averaged quantities in periodic and chaotic dynamical systems. The 
new method overcomes the difficulties that cause fundamental failure for traditional 
sensitivity analysis methods. The method is demonstrated on the periodic van der Pol 
oscillator, the chaotic Lorenz attractor, and a simple aero-elastic oscillation model that 
exhibits mixed periodic and chaotic behavior. These applications demonstrate the 
effectiveness of our new sensitivity computation algorithm in many complex nonlinear 
dynamics regimes. These include fixed points, limit cycles, quasi-hyperbolic and non- 
hyperbolic strange attractors. 
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Appendix A. Limit of derivative versus derivative of limit. 

For a given function f[x; s) of x parameterized by s, the limit of its derivative 

lim (A.l) 

s->-oo dx 

is not always equal to the derivative of the limit 

(lim (A.2) 

An intuitive counter example relevant to this paper is the following parameterized 
function 

/(a;,s) = ^sine^^a; (A.3) 



Because the derivative of the function 

df 



, (x, s) = — cos e^'^x (A.4) 
dx s 
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goes to infinity as s — >■ cx), the limit of its derivative 

lim ^ (A.5) 

does not exist. However, tlie limit of the function 

lim f{x,s) = ; (A.6) 

therefore, the derivative of the Hmit 

±(lijn fix,s))^0. (A.7) 
ax \s->-oo / 

Thus the limit of derivative is not equal to the derivative of limit. 

Appendix B. Least squares solution of a linear differential equation. 

The least squares solution of a linear differential equation 

min ||m||2 , s.t. —=Au + b{t), 0<t<T, (B.l) 

u{t) , — oo<t<oo clt 

automatically forces the divergent and convergent components of the solution propa- 
gate in different time directions. To illustrate this property, consider the case when 
M is a scalar. In this case, the matrix A is also a scalar A — X. Because Equation 
(jB.ip is a quadratic minimization problem with linear constraints, the principle of 
superposition applies, and it is sufficient to consider the right hand side b{t) = 6(t) 
being a Dirac delta function. The least squares solution satisfies 

du 

min ||u||2 , s.t. — = Xu + 6{t) , 0<t<T, (B.2) 

u{t), — oo<t<oo (it 



One can verify that for A > 0, the solution to Equation (IB.2I) is 



<t)=r'' (B.3) 

while for A < 0, the solution to Equation (IB.2p is 

u{t) = <^ B.4 

These solutions indicates that the impulse response of an "unstable" equation (A > 0) 
propagates backwards in time; while the impulse response of a "stable" equation 
(A < 0) propagates forward in time. 

For vector equations containing both stable and unstable modes, the least squares 
solution automatically separates each mode, allowing them to propagate in different 
time directions. We consider two examples. For simplicity consider A > 0. 



II II d,u 
mm ||u||2 , s.t. — 

u{t), — OQ<t<00 dl 



A 

-A 



5{t) 

m 



0<t<T , (B.5) 



mm \\u\\ 



S.t. 



du 




A 



A 




25{t) 




Q<t<T 



(B.6) 
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In the first example, it is easy to verify that the solution is a superposition of the 
solution of the stable and unstable scalar problem: 



I e^* i < 
[0 t>0 

'e-^* t > 
i < 



ui{t) = 

U2{t) = 

In the second example, the solution can be verified as 

Je-^* t > 

j-e-^' t>0 



(B.7) 



(B.8) 



We can see that in both cases, the impulse response of the system is nonzero in both 
forward and backward time directions. Different modes corresponding to different 
rate of growth or decay propagates in different directions. 



